See Also: Tutorial page
Refer to Understanding Coordinate Systems on the Tutorial page for general information about coordinate systems.
The simplest way to describe a coordinate system is by using a coordinate system string, which in Geosoft as the following form:
"horizontal_datum / map_projection [vertical_datum]"Horizontal datum names recognized by Geosoft are listed in reference table: C:\Program Files\Geosoft\Desktop Applications 9\csv\datum.csv. The horizontal_datum, map_projection and vertical_datum should use the well-known names defined by the EPSG Geodetic Parameter Registry.
Map projections recognized by Geosoft are listed in reference table: C:\Program Files\Geosoft\Desktop Applications 9\csv\transform.csv
The vertical_datum string can currently be any descriptive string, including "geoid" or "geodetic", though to future-proof your code we recommend that you use the common short name for a known vertical datum. For example, the "North American Vertical Datum of 1988" is commonly called "NAVD88".
To create an instance of a coordinate system you can provide the coordinate system string to <geosoft.gxpy.coordinate_system.Coordinate_system()> to create a Coordinate_system instance. For example, the following code defines the "UTM zone 15N" projection on the "NAD83" datum, and the "NAD27" geographic coordinate system. These defined systems are then used to translate spatial coordinates from one system to the other.
In [1]:
import geosoft.gxpy.gx as gx
import geosoft.gxpy.coordinate_system as gxcs
import numpy as np
# create context
gxc = gx.GXpy()
# define coordinate systems and a transformer
cs_utm = gxcs.Coordinate_system('NAD83 / UTM zone 15N')
cs_nad27 = gxcs.Coordinate_system('NAD27')
cs_transform = gxcs.Coordinate_translate(cs_utm, cs_nad27)
# example transform a single (x, y) coordinate
lon_lat = cs_transform.convert((345000, 64250000))
print('(lon, lat): {}'.format(lon_lat))
# example transform a single (x, y, elevation) coordinate
print('(lon, lat, elevation): {}'.format(cs_transform.convert((345000, 64250000, 50))))
# example translate a list of (x, y, z) tuples
locations = [(345000., 64250000, 50), (345500, 64250000, 60), (346000, 64250000, 70)]
nad27_locations = cs_transform.convert(locations)
print('Geographic: {}'.format(cs_nad27))
for xyz in nad27_locations:
print(xyz)
# example transform a numpy array in-place
data = np.array([[345000, 64250000, 50, 55000],
[345500, 64250000, 60, 55150],
[346000, 64250000, 70, 56000]],
dtype=float)
nad27_locations = cs_transform.convert(data, in_place=True)
for xyz in data:
print(xyz)
# compare coordinate systems
print(cs_utm == cs_nad27)
print(gxcs.Coordinate_system('WGS 84') == gxcs.Coordinate_system('WGS 84'))
print(gxcs.Coordinate_system('GDA94 [geodetic]') == gxcs.Coordinate_system('GDA94 [geoid]'))
All spatial information in the Geosoft environment will have a defined coordinate system, which if not explicitly set or inherited will be named "*unknown". Geosoft functions that mix spatial data from different coordinate systems will automatically transform coordinates to be in the coordinate system required by the context. Spatial locations that have an "*unknown" coordinate system are generally assumed to match the context in which the coordinates are used.
When data is imported into Geosoft it is good practice to define the coordinate system my assigning the known coordinate system to the coordinate_system property of the data.. If you import data via one of the Geosoft import functions, and the data has a way of describing its coordinate system, the Geosoft import function will make a best effort to set the the coordinate_system property from the data. For example, ESRI data files often have a well defined coordinate system which is recognized when the data is imported.
The geosoft.gxpy module exposes various spatial data structures via classes, such as the Geosoft_database class, the Grid class, the Geometry class, and the View class. All spatial classes in the Geosoft environment will have a property named coordinate_system, which is an instance of the geosoft.gxpy.coordinate_system.Coordinate_system class. The coordinate_system property can be used get or set the instance coordinate systems using any form supported by the Coordinate_system class constructor. The following script demonstrated the versatility of how to establish a coordinate systems using various forms. These examples are shown applied to a geosoft.gxpy.grid.Grid instance, but may be equally applied to any spatial class instance.
In [2]:
import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrd
gxc = gx.GXpy()
# create a memory grid as an example of a spatial object
grid = gxgrd.Grid.new(properties=({'nx': 10, 'ny': 10}))
print('Initially:', grid.coordinate_system) # initially '*unknown"
# define by a Geosoft-style coordinate system name. Parameters are derived from internal Geosoft tables.
grid.coordinate_system = "NAD83 / UTM zone 17N"
print('\nFrom name:', grid.coordinate_system)
print(grid.coordinate_system.gxf)
print(grid.coordinate_system.coordinate_dict())
# example use of GXF strings to change the datum to NAD27. Here we remove the name and local datum transform
# and allow the Coordinate_system class to complete parameters for NAD27 from the tables.
gxf = grid.coordinate_system.gxf
gxf[0] = ''
gxf[1] = "NAD27"
gxf[4] = ''
grid.coordinate_system = gxf
print('\nFrom gxf:', grid.coordinate_system.gxf)
# fully explicit definition of UTM zone 17N on NAD27 datum using GXF strings.
grid.coordinate_system = ['',
'NAD27',
'"Transverse Mercator",0,-87,0.9996,500000,0',
'm,1',
'"*local_datum",-8,160,176,0,0,0,0']
print('\nFrom gxf:', grid.coordinate_system.gxf)
# ... from a json string. Note how to properly escape the string embedded in a string.
js = '{"units": "m,1", "datum": "NAD27", "projection": "\\"Transverse Mercator\\",0,-87,0.9996,500000,0"}'
grid.coordinate_system = js
print('\nFrom json:', grid.coordinate_system.gxf)
# ... from an ESRI WKT string
wkt = 'PROJCS["NAD_1927_UTM_Zone_16N",' + \
'GEOGCS["GCS_North_American_1927",' + \
'DATUM["D_North_American_1927",' + \
'SPHEROID["Clarke_1866",6378206.4,294.9786982]],' + \
'PRIMEM["Greenwich",0.0],' + \
'UNIT["Degree",0.0174532925199433]],' + \
'PROJECTION["Transverse_Mercator"],' + \
'PARAMETER["False_Easting",500000.0],' + \
'PARAMETER["False_Northing",0.0],' + \
'PARAMETER["Central_Meridian",-87.0],' + \
'PARAMETER["Scale_Factor",0.9996],' + \
'PARAMETER["Latitude_Of_Origin",0.0],' + \
'UNIT["Meter",1.0],' + \
'AUTHORITY["EPSG",26716]]'
grid.coordinate_system = wkt
print('\nFrom wkt:', grid.coordinate_system.esri_wkt)